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Power Law Distributions of Seismic Rates 



O 

o 

(N 
o 

Q 

0^ 



Ph. 

O 
D 

bp; 

O 
Oh 



> 

o 
o 

(N 

O 
O 

Oh- 



A. Saichev^'^ and D. Sornette^^^ 

^Mathematical Department, Nizhny Novgorod State University, 
Gagarin prosp. 23, Nizhny Novgorod, 603950, Russia 
^Institute of Geophysics and Planetary Physics, University of California, Los Angeles, GA 90095 
^Institute of Geophysics and Planetary Physics and Department of Earth and Space Sciences, 
University of California, Los Angeles, CA 90095 
^Lahoratoire de Physique de la Matiere Condensee, 
CNRS UMR 6622 and Umversite de Nice-Sophia Antipolis, 06108 Nice Cedex 2, Franc^ 

(Dated: February 2, 2008) 

We report an empirical determination of the probability density functions Pdata('') of the number 
r of earthquakes in finite space-time windows for the California catalog. We find a stable power law 
tail Pdata(?') ^ l/r^"*"^ with exponent « 1.6 for all space (5 x 5 to 20 x 20 km^) and time intervals 
(0.1 to 1000 days). These observations, as well as the non-universal dependence on space-time 
windows for all different space-time windows simultaneously, are explained by solving one of the 
most used reference model in seismology (ETAS), which assumes that each earthquake can trigger 
other earthquakes. The data imposes that active seismic regions are Cauchy-like fractals, whose 
exponent 5 = 0.1 ± 0.1 is well-constrained by the seismic rate data. 
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Seismicity is perhaps the best example of a self- 
organizing process exhibiting scahng diagnosed with so 
many power laws: the Gutenberg-Richter distribution 
- l/£'i+'5 (with /3 « 2/3) of earthquake energies E; the 
Omori law ^ Xjt^ (with p « 1 for large earthquakes) 
of the rate of aftershocks as a function of time t since 
a mainshock; the productivity law ~ i?" (with a k. 2/3) 
giving the number of earthquakes triggered by an event of 
energy E hi ; the power law distribution ^ 1 jL? of fault 
lengths L Q ; the fractal structure of fault networks Q 
and of the spatial organization of earthquake epicenters 
0; the distribution l/s^+'' (with 5 > 0) of seismic stress 
sources s in earthquake focal zones due to past earth- 
quakes Related universal laws for the distribution 
of waiting times and seismic rates between earthquakes 
have recently been derived from the analyses of space- 
time windows 0- 

Here, we report and explain theoretically an addition 
empirical power law: the numbers r of earthquakes in fi- 
nite space-time windows for the California SCEC catalog, 
over fixed spatial boxes 5x5 km^ to 20 x 20 km^ and time 
intervals dt — \, 10, 100 and 1000 days are distributed 
with a stable power law tail PdataC?") ~ Xjr^^^ with ex- 
ponent [I ~ 1.6 for all time intervals. This result has im- 
portant implications in constraining the physics of earth- 
quakes and in estimating the performance of forecast- 
ing models of seismicity. For the former, we show that 
this observation can be rationalized by a simple stochas- 
tic branching model (ETAS) which uses the above men- 
tioned Gutenberg-Richter, Omori, and productivity laws 
applied to a fractal spatial geometry of earthquake epi- 
centers. The fundamental physical ingredient is that each 
earthquake can trigger other earthquakes ( "aftershocks" ) 
and an earthquake sequence results in this model from 
the cascade of aftershocks triggered by each past earth- 



quake. In addition, the growing efforts in earthquake 
forecasts requires estimating the performance of a fore- 
casting model by a likelihood function, which are cur- 
rently based on Poisson probabilities calculated using 
declustered catalogues. Our work shows that sponta- 
neous fluctuations of the number of triggered earthquakes 
in space-time bins, due to the cascades of triggering pro- 
cesses, may lead to dramatic departures from the Poisson 
model used as one of the building block of standard test- 
ing procedures. 
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FIG. 1: Empirical probability density functions (pdf) 
P&!Az.(r) of the number r of earthquakes in the space-time 
bins of size 5x5 km^ and dt = \ days of the SCEC cata- 
log. The straight line is the best fit with a pure power law 
Q. The continuous line is the theoretical prediction derived 
below for dt — 1 days and with no adjustable parameters, as 
the parameters are fixed to the values adjusted from the fit 
for dt = 1000 days for which the pdf is strongly curved. 
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In order to maximize the size and quality of the data 
used for our analysis, we consider the time interval 
1994 - 2003 in a region from approximately 32° to 37°N 
in latitude, and from —114° to —122° in longitude, of 
the Southern Californian earthquakes catalog with re- 
vised magnitudes Ml > 1-5, which contains a total of 
86,228 earthquakes. The spatial domain is covered by 
square boxes of (L = 5 km) x (L = 5 km). Other larger 
box sizes given similar results. We consider time win- 
dows from dt = 1 day to dt = 1000 days. Figure ^ plots 
the empirical probability density functions Pdatair) of the 
number r of earthquakes in the space-time bins described 
above for dt = 1 days. The straight line is the best fit 
with a pure power law 

Pdata(r) ^ l/ri+'^ (1) 

over the range 1 < r < 100. Similar power law tails 
are present in the tail for the other time windows. The 
estimation for fi is stable since the fitted values are fi ~ 
1.65 for dt — 100 days, ^ = 1.75 for dt = 10 days and 
/i = 1.60 for dt — 1 days. However, the pdf becomes 
more and more curved in the larger portion of the bulk 
as the size dt of the time window is increased 0. This 
behavior can be explained by the theory described below. 

The ETAS (Epidemic- Type Aftershock Sequence) 
model of triggered seismicity is based on the three first 
well-founded empirical laws mentioned above. Its ap- 
peal lies in its simplicity, its power of explanation of 
other empirical observations (see for instance 0| and ref- 
erences therein) and its wide use as a benchmark. The 
ETAS model belongs to the general class of branching 
processes with infinite variance of the number of proge- 
nies per mother, with a long-time (power law) memory of 
the impact of a mother on her first-generation daughters 
described by the empirical Omori law for aftershocks. 
These two ingredients together with the mechanism of 
cascades of branching have been shown to give rise to 
subdiffusion and to non mean-field behavior in the distri- 
bution of the total number of aftershocks per mainshock 
, in the distribution of the total number of generations 
before extinctions and in the distribution of the total du- 
ration of an aftershock sequence before extinction |0| . 

In the ETAS model, each earthquake is a potential pro- 
genitor or mother, characterized by its conditional aver- 
age number Nm = KiJ.{m) of children (triggered events or 
aftershocks of first generation), where /i(m) — 10"(™~'"o) 
is proportional to the average productivity of an earth- 
quake of magnitude m > mo k is a constant fac- 
tor and niQ is the minimum magnitude of earthquakes 
capable of triggering other earthquakes. For a given 
earthquake of magnitude m and therefore of mark ii{m), 
the number r of its daughters of first generation are 
drawn at random according to the Poissonian statistics 
Pn{r) = e^^'" = -^^^^r- e"''^. The challenge of our 
present analysis is to understand how the exponential 
Poisson distribution is renormalized into the power law 



by taking into account all earthquake triggering paths 
simultaneously over all possible generations. The ETAS 
model is complemented by the normalized Gutenberg- 
Richter (GR) density distribution of earthquake magni- 
tudes p(m) = b In(lO) 10-''(™-'"°),m > mo. This mag- 
nitude distribution }3(m) is assumed to be independent of 
the magnitude of the triggering earthquake, i.e., a large 
earthquake can be triggered by a smaller one. Combin- 
ing the GR and the productivity laws shows that the 
earthquake marks /i and therefore the conditional aver- 
age number Nm of daughters of first generation are dis- 
tributed according to the normalized power law 

Pp(m) = nrpr , 1 < A* < +00, ■y = b/a. (2) 

For earthquakes, 6 « 1 and 0.5 < a < 1 giving 1 < 7 < 2. 
This range 1 < 7 < 2 implies that the mathematical 
expectation of /x and therefore of TV™ (performed over 
all possible magnitudes) is finite but its variance is infi- 
nite. Given 7, the coefficient k then controls the value of 
the average number n (or branching ratio) of children of 
first generation per mother: n = (Nm) = ^(/i) = n-^, 
where the average (N^) is taken over all mothers' mag- 
nitudes drawn from the GR law. Recall that the val- 
ues n < 1, n = 1 and n > 1 correspond respectively 
to the sub-critical, critical and super-critical branching 
regimes. The last ingredient of the ETAS model con- 
sists in the specification of the space-time rate func- 
tion Nm ^{r — fijt — ti) giving the average rate of first 
generation daughters at time t and position r created 
by a mother of magnitude m > mo occurring at time 
ti and position r^. We use the standard factorization 
^{x,t) = ^{t)(f){x). The time propagator <i>(t) has the 
Omori law form $(<) = (c+ty+o ^(^) where H{t) is 
the Heaviside function, < 9 < l,cisa regulariz- 
ing time scale that ensures that the seismicity rate re- 
mains finite close to the mainshock. The space propa- 
gator is 4){x) = 2^(a:^-Hd^)'('i+2)/2 ■ Thc ucxt ingredient of 
the ETAS model is to assume that plate tectonic motion 
induces spontaneous mother earthquakes, which are not 
triggered by previous earthquakes, according to a Pois- 
sonian point process, such that the average number of 
spontaneous mother earthquakes per unit time and per 
unit surface is g. In the ETAS branching model, each 
such spontaneous mother earthquake then triggers inde- 
pendently its own space-time aftershocks branching pro- 
cess. The last ingredient of our theory is to recognize 
that, at large scale, earthquakes are preferentially clus- 
tered near the plate boundaries while, at smaller scales, 
earthquakes are found mostly along faults and close to 
nodes between several faults. We thus extend slightly 
the ETAS model to allow for the heterogeneity of the 
spontaneous earthquake sources g reflecting the influence 
of pre-existing fault structures, some rheological hetero- 
geneity and complex spatial stress distributions. For this, 
we use the distribution of the stress field in heteroge- 
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neous media and due to earthquakes j5( which is found 
close to a Cauchy distribution. The simplest prescription 
is to assume that g is itself random and distributed ac- 
cording to / (^■(fy^ 7 where {g) is then statistical aver- 
age of the random space-time Poissonian source intensity 
g. In the numerical applications, we shall use the form 
/,(x) = ^(l + f)-'-' (<5>0). 

Due to the independence between each sequence trig- 
gered by each spontaneous event, the generating proba- 
bility function (GPF) of the number of events (including 
mother earthquakes and all their aftershocks of all gener- 
ations), falling into the space-time window {[t, i-f t] x S} 
is equal to 

e,p(0,r,5) =e-^^(^'^"5) (3) 

OO 

where L{z,t,S) = dt JJ dx[l ~ 9(z, t, r, 5; a^)] + 

— oo 

oo 

/([di // dx[l-e{z,t,S;x)][l-Isix)]+ dt JJ dx[l- 

-oo S 

zQ(z, t, iS; x)]. The first summand in L(z, r, S) describes 
the contribution to the GPF Qsp from aftershocks trig- 
gered by mother earthquakes that occurred before the 
time window (i.e. at instants t' such that t' < t). The 
corresponding GPF Q{z,t — t',T,S;x) of the number 
of aftershocks triggered inside the space-time window 
{[t,t + t] X S} by some mother event that occurred at 
time t' — satisfies the relation 

e(z, t, r, S; x) = G[l - ^{z, t, r, S; x)] , (4) 

where the auxiliary function ^!{z,t,T,S;x), describing 
the space-time dissemination of aftershocks triggering by 
some mother event, is equal to 

*(z, t, r, 5; x) = ^{x, i) [1 - e(z, t, t, 5; a;)] 

+<P{x,t + T)(g)[l-e{z,T,S;x)]+ (5) 
(1 - z)<P{x, t + r) ® Is{x)Q{z, T, S; x) . 

The function Is{x) in jSJ is the indicator of the space 
window S and G{z) in Q is the GPF of the num- 
ber i?i given by 10) of first generation aftershocks trig- 
gered by some mother earthquake, given by G{z) — 
jK^{l — z)'^ r(— 7, k(1 — z)). The last two summands of 
L(z,T,S) describe the contribution of aftershocks trig- 
gered by earthquakes, occurring inside the time window 
(i.e., t' G [t,t + T]). The second (resp. third) term cor- 
responds to the subset spatially outside (resp. inside) 
the domain S. These last two terms depend on the GPF 
<d{z,T,S;x) = <d{z,t = 0,t,S;x) of the numbers of af- 
tershocks triggered till time r inside the space window S 
by some mother event arising at the instant t — and at 
the point x. It follows from |0J and Q that it satisfies 
the relations 

e(z,r,5;a;) = G[l - *(z,r,5;a;)] (6) 



and 

*(z, T, S; x) = $(a;, r) O [1 - e(z, r, 5; x)] 

+ (1 - z)^{x, t) ® Is{x)e{z, T, 5; x). ^'^ 

In addition, we shall need the GPF 

Q{z,S;x) = Q{z,T = oo,S;x) (8) 

of the total numbers of aftershocks triggered by some 
mother earthquake inside the area S. As seen from © 
and (TJ, it satisfies the relations 

e(z, 5; x) ^ G[l - *(z, S; x)] (9) 

and ^i{z, S;x) = I - 4>{x) e(z, S; a;) + (1 - z)(j){x) (g) 
Is{x)Q{z,S;x). 

Taking into account the distribution of the source in- 
tensities g amounts to averaging equation (|3Jl over g 

weighted with the statistics '{^ / ("(fy) • This gives 

e,p(z,T;5) = /[(e)L(z,T,5)], (10) 

where f{u) is the Laplace transform of the pdf f{x). 

To go further, we make two approximations. If the 
time duration r of the space-time window is sufficiently 
large, it can be shown that the statistical averages of 
the seismic rates become independent of r. It seems 
reasonable to conjecture that the GPF Q(z,t,S;x) of 
the total number of aftershocks triggered by some earth- 
quake source inside the space domain S until time r coin- 
cides approximately with the saturated GPF Q{z, S; x) of 
the total number of aftershocks triggered by some earth- 
quake source inside the space domain S. Within this 
approximation of large time windows, the effect of af- 
tershocks triggered by earthquake sources occurring till 
the beginning t of the time window is negligible. Within 
this large time window approximation, one may ignore 
the first term in the contribution to Q and replace 
Q{z^t,S;x) by Q{z,S;x) in the remaining terms. As a 
result, L(z,T,S) in Q takes the following approximate 

oo 

form L{z,T,S) ~ t//[1 — Q{z,S;x)][l — Is{x)]dx + 

— oo 

T JJ[1 — zQ{z,S;x)]dx, where Q{z,S;x) is the solution 
s 

of e = G [6 (g) (/) - (1 - 2)756 ® 0] . To find a reasonable 
approximate expression for the sought GPF Q{z,S;x), 
notice that if the spatial extend £ of the window is larger 
that the characteristic scale d of the space kernel, or if 
n is close to 1, then the characteristic spatial scale as- 
sociated with the GPF Q{z,S;x) becomes greater than 
d. Therefore, without essential error, one may replace 
Q (i) (phy &. In addition, we take into account the finite- 
ness of the domain S by using the factorization proce- 
dure: Is{x)<d{z,S;x) (j){x) ~ Q{z,S;x) ps{x), where 
Ps{x) remains to be specified. This amounts to replacing 
a convolution integral by an algebraic term. This factor- 
ization approximation is a crucial step of our analysis 



4 



and is justified elsewhere 0|- As a result of its use, the 
nonlinear integral equation for &{z, S; x) transforms into 
the functional equation 6 = G[{l + {z — l)ps{x))Q]. This 
approximation leads to (R) = j^PSj where (i?) is the 
average of the total number of events in the space-time 
window. The effective parameter can be determined 
from the consistency condition such that (R) be equal to 
the true {R{S; x)), which can be calculated exactly. This 
gives 



T) ~ ~ 1 — ?^ 

-{R{S;x)), ps{q) ^ Is{q)HQ)- 



n<f){q) 
(11) 

For £ ^ d, the factor ps{x) approaches a rectangular 
function. We can use this observation to help determine 
the statistics of the number of events in a finite space- 
time window, using the approximation (a;) ~ const ~ 
p for X ^ S. We define the parameter p as the space 
average of ps{x) over the window's area S: 



This approximation allows us to get 

T JJ [1 - ze{z,S;x)]dx ~ tS[1 - ze{z;p)] 
s 



(12) 



(13) 



where Q{z,p) is the solution of Q{z;p) = G[(l + (z — 
1)pMz;p)]. 

Complementarily, a study of ps{x) shows that it is 
small outside the window space domain S. This implies 
that, outside 5, one may replace the functional equation 
on 6 by 1 — 9(z, S; x) ~ jT^II — z)ps{x). Therefore, we 

oo 

get r // [1 - e(z,5; x)][l - Is{x)]dx ^ qrS ^(1 - z), 

— oo 

oo 

where g = ^ JJ ps{x)[l — Is{x)]dx. Taking into account 

— oo 

oo 

that // ps{x)dx = S, we obtain q ~ 1 — p. Putting all 

— oo 

these approximations together allows us to rewrite the 
expression of L{z,t,S) in Q as 



-{l-p)(l- z) + l- zQ{z;p) 



L{z, T, S) ~ tS 

(14) 

The factorization procedure obtains the characteristic 
features of the space-time branching process in a finite 
space-time domain, at the cost of an adjustable parame- 
ter p. 

Starting from the general expression H10|l of the GPF 
0sp(z, t;<S) with the approximation (|14|l for L[z,t,S) 
and using the relationship between the probability 
Psp{r;p,p) and its GPF in the form of its integral rep- 
resentation, we obtain the following expression valid in 
the limit of sufficiently large time windows Pspir; p,p) = 



^JHp [t^(1 -p)(1 - ^) + 1 - ^e(z;p)l) ^. 
Introducing the new integration variable y = (1 + (z — 
l)p)e(z;p) ^ ^^z(y)-i(^+p-l),by 
construction of y, 6(2:; p) = G(y) which allows us to ob- 
tain the following explicit expression Psp(^; PiV) — 2^ x 

dZ(y) dy_ 



§nP T^(i-p)(i-^(?y)) + i-^(y)G(y) 



This expression makes a precise quantitative prediction 
for the dependence of the distribution Psp(r;p,p) of 
the number r of earthquakes per space-time window as 
a function of r, once the following model parameters 
are given: the branching ratio n, the exponent 7 of 
the distribution of productivities, the exponent 5 of 
the distribution of spontaneous earthquake sources, the 
fraction ps of direct (first generation) aftershocks which 
fall within the domain <S, and the average number p 
of spontaneous earthquake source per space-time bin 
defined hy p = (g) tS. The theoretical curve in Fig. Q] 
is obtained by a numerical integration of Psp{r; p,p) = 
for the set of parameters n — 0.96, 7 = 1.1, ps — 0.25, 
6 = 0.15 and p = 0.0019 dt with dt in units of days 
(thus equal to 100 for Fig. H)- These parameters give 
the best fit for the large time window dt = 1000 days. 
They have been kept fixed for the other time windows 
which exhibit very different shapes in their bulk. The 
theory is thus able to account simultaneously for all the 
considered time windows, with no adjustable parameters 
for the three smallest time windows |2|. 

We would also like to stress that, according to our the- 
ory, the value of the exponent /i w 1.6 used in to 
fit the tails of the distributions of seismic rates is de- 
scribing a cross-over rather than a genuine asymptotic 
tail. Recall that the distribution of the total number 
of aftershocks has two power law regimes ^ for 
r < r* ~ 1/(1 - n)'^/(''~i) and - l/r^+'' for r > r* 0]. 
The existence of this cross-over together with the con- 
cave shape of the distribution at small and intermediate 
values of r combine to create an effective power law with 
an apparent exponent p ~ 1.6 larger than the largest 
asymptotic exponent 7. We have verified this to be the 
case in synthetically generated distributions with gen- 
uine asymptotics exponent 7 = 1.25 for instance, which 
could be well fitted by « 1.6 over several decades. We 
note also that Pisarenko and Golubeva with a dif- 
ferent approach applied to much larger spatial box sizes 
in California, Japan and Pamir-Tien Shan, have reported 
an exponent p < 1 which could be associated with the 
intermediate asymptotics characterized by the exponent 
1/7 < 1, found in our previous analysis By using 

data collapse with varying spatial box sizes on a Califor- 
nia catalog. Corral finds that the distribution of seismic 
rates exhibits a double power-law behavior with p ^ Q 
for small rates and p 1.2 for large rates 0]. The first 
regime might be associated with the non universal bulk 
part of the distribution found in our analysis. The second 
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regime is compatible with the prediction for the asymp- 
totic exponent = 7. In conclusion, we have offered a 
simple explanation of the power law distribution of seis- 
mic rates, which is derived from the other known power 
laws and the physics of cascades of earthquake triggering. 
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